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Abstract 

Deriving the Feynman rules for lattice perturbation theory from actions and op- 
erators is complicated, especially when improvement terms are present. This phys- 
ically important task is, however, suitable for automation. We describe a flexible 
algorithm for generating Feynman rules for a wide range of lattice field theories in- 
cluding gluons, relativistic fermions and heavy quarks. We also present an efficient 
implementation of this in a freely available, multi-platform programming language 
(Python), optimised to deal with a wide class of lattice field theories. 
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1 Introduction 



Non-abelian quantum field theories such as QCD are believed to explain much 
of particle physics, at least at energy scales probed by current particle accel- 
erators. Perturbative expansions of the theory do not, however, converge at 
hadronic energy scales. That, and the belief that non-perturbative physics may 
also contribute to certain states, makes the lattice regularisation of quantum 
field theories extremely important. Inherently non-perturbative calculations 
can then be carried out using Monte Carlo simulations. 
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Dividing space and time into a grid with lattice spacing a, however, excludes 
ultraviolet modes with momenta of 71/ a or higher. A renormalisation pro- 
gramme is therefore required to connect lattice measurements to their con- 
tinuum counterparts. Such renormalisation factors are particularly important 
for QCD matrix elements and fixing the couplings and masses present in the 
Lagrangian. Renormalisation is also needed to determine the strong coupling 
as and to relate the lattice regularisation scale Aiat to the more familiar Aqcd- 
It is also used to "improve" the lattice actions in an attempt to reduce the 
discretisation errors at given lattice spacing. 

In a limited number of cases the renormalisation constants can be determined 
using non-perturbative techniques. Results at finite lattice spacing, however, 
can depend upon the method used (e.g. [1]), and non-perturbative methods 
do not cope well with mixing of operators under renormalisation. For these 
reasons there is a strong interest in lattice perturbation theory. 

Given that perturbation theory fails in low energy QCD, we may ask why 
it should work on the lattice. An argument for its use is given in [2]: the 
renormalisation factors may be thought of as compensating for the ultraviolet 
modes excluded by the lattice regulator. For typical lattices a < 0.1 fm, and 
the excluded modes have momenta in excess of 5 GeV. At these scales the 
running QCD coupling as is small enough that perturbation theory should 
rapidly converge. The wide range of results recently reviewed in [3,4] show 
perturbation theory can be used for a large range of lattice QCD processes. 
It is an assumption that non-perturbative effects do not contribute on these 
short length scales. In a few cases we can test this directly by comparing high 
order perturbative calculations with Monte Carlo simulations at a range of 
weak couplings [5,6,7,8,9]). The non-perturbative contributions to the studied 
quantities are very small. Other comparisons, such as [1], cannot distinguish 
non-perturbative effects from higher loop perturbative corrections. It there- 
fore remains that lattice perturbation theory provides the only systematically 
improvable method for determining the full range of renormalisation constants 
[3]. 

As in the continuum, the calculation of lattice Feynman diagrams is a two 
stage process. The lattice action and operators must first be Taylor expanded 
to give the propagators and vertices that form the Feynman rules (which we 
refer to as the "vertex expansion" stage). Following this, these rules must be 
used to construct and evaluate Feynman diagrams, possibly after algebraic 
simplification (the "Feynman diagram evaluation" stage). 

The main obstacles in the latter task are the presence of Lorentz symmetry 
violating terms at finite lattice spacing and the complications of replacing 
momentum integrals by discrete sums. The calculations are therefore usually 
done using computer programs hke Vegas [10], Form [11] or other propri- 
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etary mathematical packages. 

Expanding the lattice action and operators to obtain Feynman rules is far more 
complicated than in the continuum. Firstly, lattice gauge fields are elements of 
the Lie group rather than the algebra of the gauge group. We must therefore 
expand exponentials of non-commuting fields to obtain the Feynman rules. 
Secondly, modern lattice theories contain a large number of irrelevant (in the 
renormalisation group sense of the word) terms chosen to improve specific 
aspects of the Monte Carlo simulation, such as the rate of approach to the 
continuum or chiral limits of QCD. 

There is, however, no unique prescription for these terms, and the choice de- 
pends on that quantities we are most interested in simulating. As a result, a 
large number of actions and operators are currently in use. Although the dif- 
ferences may be subtle, each choice provides a separate rcgularisation of QCD 
with its own set of renormalisation constants and, most relevantly here, Feyn- 
man rules. At present the complications of the expansions has meant that the 
availability of renormalisation factors has lagged far behind developments in 
lattice improvement. In many cases this has restricted the physical predictions 
obtained from the simulations. 

As a result, there is a strong need for an automated method for deriving 
lattice Feynman rules in a fiexible way for a range of different theories. The 
generation should be rapid enough not to constrain our choice of action, and to 
avoid errors we should be able to specify the action in a compact and intuitive 
manner (such as using nested link smearing prescriptions). The evaluation of 
the Feynman diagrams can be computationally intensive, and may be carried 
out on costly supercomputing facilities. Parsimony and software availability 
dictate that the rules should be separately calculable in advance, and rendered 
in a machine readable format that can be copied to any computer for later 
Feynman diagram evaluation. 

In this paper we describe such a method. 

Automated expansion of lattice actions^ is not a new concept, having been 
described for gluonic actions by Liischer and Weisz in 1986 [12]. An imple- 
mentation of this has been used in [13,14,15]. A similar method is employed in 
[16]. We present here a new algorithm suited to expansion of not only gluonic 
actions, but also those of complicated relativistic fermionic actions and heavy 
quarks, such as in NRQCD. As in [12], the expansion is independent of the 
boundary conditions allowing, for instance, the use of twisted boundary con- 
ditions to regulate infrared divergences in a gauge-invariant manner [17,18] 
or otherwise change the discrete momentum spectrum [19]. We also describe 

^ We shall understand the term "actions" to include measurement operators from 
now on. 
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details of an implementation of this algorithm which wc have used for calcu- 
lations of the renormalised anisotropy in gauge theories [20,21], to study the 
mean link in Landau gauge for tadpole improvement [9] and to measure the 
electromagnetic decays of heavy quark systems using NRQCD [22,23]. The 
code is flexible and can be easily extended to cope with a full range of prob- 
lems, some of which we discuss in Section 5. We are happy to share this code 
with interested readers. 

The structure of the paper is as follows. Our method clearly stems from [12] 
and in Section 2 we review their theory and notation. The algorithm itself dif- 
fers markedly from [12], not least in being able to deal with fermionic actions, 
and is described in Section 3. In Section 4 we turn to implementation of the 
algorithm, explaining the steps taken to ensure the code can cope with the 
more complicated theories. Whilst the notation is tilted towards our version in 
the Python programming language, the optimisations are clearly applicable 
to any realisation of the algorithm. We present our conclusions in Section 5. 
Technical details of the the data structures employed are relegated to Appen- 
dices. 



2 The lattice 



A cubical space-time lattice A in D dimensions consists of sites labelled by a 
vector 33 G A with components that are integer multiples of a lattice spacing a, 
which we will set to be one (a (bare) lattice anisotropy can be introduced 
through rescahng of couphng constants in the action [20]). The directions 
of the lattice axes are labelled e {1, 2, . . . , £)}. If is a right-handed 
basis set consisting of unit vectors, we define corresponding backward vectors: 



A path consisting of I links starting at site x can be specified on the lattice 
by an ordered set of signed integers, Sj e [—D, . . . , —1, 1, . . . , 

C{x, y; s) = {x, y,s= [sq, si, . . . , s/_i]} . (1) 

The j point on the path is 

{X , 7 = 0, 

(2) 
Zj^i + aes._, , 0<j<l, 

and the endpoint of the path is y = Zi. 

For a periodic lattice with sites in the /i direction (and volume V — H^u ^//) 
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the momentum vectors are 

= — f^,---,^) , 0<^^<L^, ^^ez, (3) 

a \Li Ld J 

and X^fe stands for sums over the integers The Fourier expansion of a field 

(j) is 

m = E e-^'^'Xa^) , = ^ E ■ (4) 

X k 

Different boundary conditions (e.g. twisted [12,20,9]) change the colour fac- 
tors and momentum spectrum. Since neither are used explicitly in the vertex 
expansion below, the same reduced vertex function output can be used in each 
case. 



2.1 Matter fields 



We now turn to the description of lattice fields. The notation follows [12]. 

The gauge field associated with a link is U^y{){x) G SU{N). Let U denote the 
full configuration of such links. The perturbative gauge potential associated 
with the link is defined through 



where g is the bare coupling constant. The potential G alg(S'?7(A^)) is 
associated with the midpoint of the link. Expanding in the anti-Hermitian 
generators of SU{N): 

A^ = Al T,, [T„, n] = -fabcTc, Tr (T^T,) = -| ■ (6) 

We define U_^{x) = Uj^{x — ae^). 

Quark fermion fields ^(a;) transform according to the representation chosen 
for the generators T^. From now on we assume this to be the fundamental rep- 
resentation (other choices will affect the colour factors, but not the underlying 
expansion algorithm). 

The Wilson line L{x, y, U) on the lattice associated with the path C{x, y; s) 
is a product of links 



L{x, y,U) = C:U = H {Zi) = H ^xp sgn(si)a^A|,.| ( Zi + -e^, ) . (7) 

-■-1 1=0'- \ z , 
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As all actions and operators can be written as a sum of Wilson lines (possibly 
terminated by fermion fields that are not themselves expanded), our goal is 
to efficiently render L as a Taylor series in the gauge potential in momentum 
space: 

L(a;,2/;A) = ^M: ^ ... Al\{k,) . . . AUK) x 

Vr{ki,Hi,ai;...;kr,Hr,ar)- (8) 

We can write the vertex functions Vr as 

Vr{ki, III, oi; . . . ; K, Hr, ar) = Cr{ai, ...,ar) Yf-{ki, //i; . . . ; K, Hr) (9) 
The matrix colour factor plays the role of the Clebsch-Gordan factor: 

a(ai,...,a,) = nT„, . (10) 

Up to differences in the colour trace structure of the the action (e.g. a mixed 
fundamental/adjoint gauge action, and discussed in Appendix B), the C^- are 
path independent. We can therefore represent the vertex functions more ef- 
ficiently by calculating just the expansion of the reduced vertex functions, 
Yf' (with an appropriate description of the colour trace structure where am- 
biguous) . The reduced vertex function can be written as a sum of monomials 

y/(fci, //i; . . . ; Hr) = ^ /„ exp - (fci •< + ... + fc^ • <) (11) 

n=l ^ 

For each combination of r Lorentz indices we have rir terms, each with an 
amplitude / and the locations v of the r factors of the gauge potential. To 
simplify this expression we have suppressed the dependence of /, f" and rij. 
on the Lorentz structure. To construct Y for given momenta, we apply the fe's 
to the position vectors of all monomials with the correct Lorentz indices. 

The v's have been drawn from the locations of the midpoints of the links in the 
path £. To avoid floating point ambiguities, it is therefore more convenient to 
express the components of all D-vectors as integer multiples of | (accounting 
for the factor of | in the exponent) . 

2.2 Realistic actions: the fermion sector 

We begin our discussion of realistic lattice actions with the fermion sector. 
The most general gauge- and translation-invariant action can be written as 

Sf{^, t/) = E E hyv^{x)r^W{x, y, U)^{y) . (12) 
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It consists of Wilson lines W defined by open paths W{x,y; s). Associated 
with each path is a coupling constant hy\> and a spin matrix r>v (which might 
be unity). 

Using the convention that all momenta flow into the vertex, the perturbative 
expansion is 



^^{pWfAp^^'^ 9,c; fci,/xi,ai; fc^, /x^, a^)?/''=(q) . 

p,q,h,c 

(13) 

The Euclidean Feynman rule for the r-point gluon-fermion-anti-fermion ver- 
tex is —g'^VF,r-i where the symmetrised vertex is: 

^F,r-(P, q, c; fci, jii, oi; . . . ; fc^, jir, Or) = 

^ ^ (7 • CF,r{h, c;ai,..., a^) a ■ YF,r(p, Q] fei, l^i] ■■■] K, A«r) ■ (14) 

a is an element of the permutation group of r objects, Sr, applied to the 
gluonic variables and normahsed by the factor of (r!). The reduced vertex 
Ypr — J2w hw^F]r Contributions from paths W. 

For all simple cases the Clebsch-Gordan colour factor is the matrix element: 

Cp^rib, C; Oi, . . . , Qr) = (T^^ . . . Tajbc ■ (15) 

The symmetrisation and calculation of colour factors will be carried out sep- 
arately when the vertex functions are reconstructed in a Feynman diagram 
calculation. 

The reduced vertex function has the structure: 



^F,r-(P, q, ki, //i; . . . ; K, Hr) ^ Y ^nfri X 

n=l 

expQ(p-a; + q-t/ + A;i-<H-... + /c, -O^ . (16) 

As we do not use explicit representations of the spin matrices, it is important 
that each monomial retains the correct spin dependence r„. 
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2.3 Realistic actions: the gluon sector 

A general gluonic action is 

5(V^,t/)=^^cpReTr[P(a;,a;,[/)] , (17) 

X ■p 

built of Wilson loops P defined by closed paths V{x,x;s), each with coupling 
constant c-p. The perturbative action is 

M^) = E^ E ••• E A2{k,)...Al{K)x 

' fcl,/il,ai hr i^lr -jdr 

VG,r(fel, /il, ai, . . .] kr, Air, Or) • (18) 

The Euclidean Feynman rule for the r-point gluon vertex function is {—g'^VG,r), 
and the vertex is [12] 

VG,r(fcl, yWl, Qi; • • • ; fcr, /^r, ^r) = 

— ^ (7 • CG,r(ai, . . . , a^) (7 • YG,r{ku Hi] kr, Hr) , (19) 

The reduced vertex = X^-p c-pYq^ is the sum of contributions from paths V. 
As before, the (r!) factor normalises the symmetrisation. Yq ,^ can be expanded 

as 

y|,(fci,//i;...;fc„/.r) = ^/„ exp -(fci-< + ... + fc,-0 . (20) 

n=l ^ ^ ^ 

In most cases we expect the lattice action to be real. For every monomial in 
Eq. (20), then, there must be a corresponding term 

(-l)7:exp-(^Efc^-<) ■ (21) 

We can therefore speed up the evaluation of the Feynman rules by removing 
the latter term, and replacing the exponentiation in Eq. (20) with "cos" for r 
even, and with "i sin" for r odd. Clearly we must identify to which terms this 
has been applied. This can cither be done by recognising conjugate contours 
in the action (e.g. S = | Tr[P + P^) and expanding only one, or by attaching 
a flag to each monomial to signal the reduction (as discussed in Section 4). 

If, in addition to the reality, the action has the form Eq. (17) with a single 
trace in the fundamental representation, the colour factors are 

CG,.(ai, . . . , a,) = i [Tr (T,, . . . T„J + (-1)'' Tr (T,, . . . T„J] . (22) 
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When symmetrising, a lot of the terms have similar Clebsch-Gordan factors: 

{1 for cr a cyclic permutation, 
(—1)'' for a the inversion. 

(23) 

We can therefore partly symmetrise the vertex over (the subgroup of cyclic 
permutations and inversion) at the expansion stage. The Xr(c") go into the 
amplitudes of the new terms coming from the partial symmetrisation: 

VG,r(fci,/^i,ai; fer,/^r,ar) = • C'G,r(«l, o^) x 

y^r- E cvXr{o) a -Yl, . (24) 

V 

The advantage of doing this is that many of the extra monomials are equiva- 
lent, and we can therefore cut down significantly the number of exponentiation 
operations required to construct Vcr- The number of remaining symmetrisa- 
tion steps (to be carried out in the Feynman diagram code) is the number of 
cosets in S^jZ^ (one for r < 3, three for r = 4 etc.). 

2.]^. Diagram differentiation 

There are many cases where Feynman diagrams need to be differentiated with 
respect to one or more momenta. Whilst this can be done numerically using an 
appropriately local difference operator, this can lead to numerical instabilities. 

It is clear from Eq. (11), that we can easily construct the differentiated Feyn- 
man vertex. Let the momentum component we wish to differentiate with re- 
spect to be q^. We first construct a rank r object r = [ri,...,rj.] which 
represents the proportion of momentum q in each leg of the Feynman dia- 
gram. Momentum conservation dictates Z^i^i = 0. For instance, for a gluon 
3-point function with incoming momenta (p, — p -|- 2g, — 2g), we would have 
r — [0, 2, —2]. The differentiated vertex is 

-^y/(fci, //i; . . . ; l^r)=f2 ^ (jivi.^ + ... + T^<.^) X 

exp ^ (fci • + . . . + fe, • <) (25) 
and so on for higher derivatives. We may therefore simultaneously calculate as 
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many differentials as we need for the cost of one exponentiation. If this momen- 
tnm expansion is placed into an appropriate data structure with overloaded 
operations, it is easy to create the Taylor series for a Feynman diagram by 
multiplying the vertex factors together. For examples of such codes, see [24]. 

It may also be necessary to differentiate diagrams with respect to parameters 
in the action, which may be present in different multiples in the amplitude of 
the monomials. Depending on the situation, we can separately expand parts 
of the action containing different powers of the parameter. Alternatively, we 
can use a single expansion and append a label to each monomial that records 
how many powers of the given parameter are contained in the amplitude. This 
label is then used in the parameter differentiation in the Feynman diagram 
code. 



2.5 Recursive path definitions 

So far we have assumed that the paths in the action are constructed from 
single links. This is, of course, always true but it is often more compact to 
specify the action as built from composite objects, such as smeared links, 
giving a separate prescription for the link smearing. For instance, the smeared 
link might be defined as the gauge covariant, weighted sum of the link and its 
adjoining staples: 

W^{x) = coU^{x) + C1YI U,{x)U^{x + e,)Ul{x + e^) , (26) 

where the coefficients Co,i define the smearing method and are chosen to op- 
timise certain aspects of the Monte Carlo simulation. The smearing may also 
be defined recursively, with smeared links inserted into additional smearing 
recipes. Examples include the "HISQ" improved staggered fermions discussed 
in [25], where the hnks are first "FAT7", then "ASQ" smeared. Rather than 
multiplying out the paths in the action to give a large sum of tangled paths 
built from links, it is more convenient to reflect the nested improvement struc- 
ture in the expansion algorithm itself. 

This is done by first defining the mapping U ^ 1^ as a sum of paths as in 
Eq. (1). For instance, Eq. (26) is represented by path {a;, a; + e^; [/x]} with 
coupling Co plus {tc, £C + e^; [z/, yU, — z/]}, {a;, £C + e^; [— z/, /x, z/]} with coupling 
Ci. Call the smearing definitions W. An action is then compactly specified by 
a path V of composite objects. These are in turn are defined by an ordered 
list [Wi, yVi, • ■ ■ , Wn] representing each step of nested improvement. The full 
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field is then defined by the recursion 



W2 : W^{U) = W2{Wi{U)) , 
W3 : W2{U) = Ws{W2{U)) , 

; (27) 



3 An expansion algorithm 



In this section we present a new algorithm for carrying out the Taylor expan- 
sion of lattice actions in a manner suited to computer implementation. 

We start by defining an object that represents a single term in the Taylor 
expansion in Eq. (11). We call this an "entity" E, and it is an ordered list: 

E = (/xi, . . . ,iir;x,y;vi,. . . , v^, /) . (28) 

The order of the entity is r. For instance, a single link comes from a path 
C — {0, e^; [//]}, and the r^^ term in its expansion in Eq. (5) is represented as 

Er = {n,...,li]0,2ei,;e^,...,e^;l). (29) 

r terms r terms 

Note that in units of | the endpoint is 2e^ and the midpoint e^. The reduced 
vertex function is a set of all the entities of all orders in the expansion, and 
we call this a "field" F = {E}. 

In practise we build a Wilson line by concatenating smaller paths (the smallest 
being the hnk). We therefore define the multiphcation of two fields so as to 

give the Taylor expansion of the resulting, longer contour. The product is 
therefore the ordered product of each entity from the first contour with every 
entity from the second: 

F{£^*£2) = F{Ci)*F{C2) , 

= {e' * E^ y E' e F{£i), E^ e F{£2)} ■ (30) 

To keep gauge covariance, entity E^ must be translated to start at the end 
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point of E^. The order of the product entity is the sum of those of the con- 
stituents. Further details are given in Appendix C. 

Addition of fields should represent the expansion of a sum of gluonic paths. 
We therefore simply combine the lists of entities: 

F{Ci + C2) = F{C,) + F{C2) . (31) 

In general, F will be a redundant representation of the polynomial containing 
two or more equivalent entities. As each entity represents a monomial which 
requires a computationally expensive exponentiation, we construct a compres- 
sion operation which compares all pairs of entities in F and combines them if 
they are equivalent: 

[E\E^] — >[E] if E' = E^ 
where E = . . . , f4 ;x',i/;vi . . . ,vl ; f + f) . (32) 

Assuming we have a translationally equivalent theory, entities need only be 
equivalent up to translation by a constant vector. For details, see Eq. (C.4). 



4 A practical implementation 

In this section we describe an implementation of the algorithm in a program- 
ming language called Python. The Python interpreter is freely available 
for a wide range of computational platforms at [26]. The complexity of some 
physical actions (notably high order NRQCD) require the implementation to 
be CPU and memory efficient. This can be achieved in PYTHON with appro- 
priate programming techniques, and without sacrificing the object orientation 
and superior list handling features of the language. 

As a first step, we choose a maximum order for the Taylor expansion. Any 
entity of higher order is discarded. The entities and all sub-lists (including 
vectors) are encoded as "tuples", which are immutable list objects to which 
the native hash function can be applied. 

To minimise the size of the dictionary, it is important that the field does not 
contain entities that are equivalent. We choose data structures for the entity 
and field specifically to prevent this. Firstly, we exploit translational invariance 
of the action to arrange that all paths start at the origin {x — 0). Entity 
equivalence in Eq. (C.4) then follows from an item by item tuple comparison. 

This comparison is most efficiently done using a hash table (i.e. an associative 
array). The Python "dictionary" is a native implementation of this, where 
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a "key" indexes an "object". In our implementation, each dictionary entry 
represents a monomial in the reduced vertex function. The key is a list of 
all information in the entity bar the amplitude, which becomes the object 
indexed by this key. Searching for equivalent entities now amounts to enquiring 
if the key already exists in the dictionary. As all key entries are integer, we 
do not need to worry about machine precision issues. We can significantly 
speed up the hashing by omitting the now redundant x = from the entity 
data structure. Independently, a moderate performance gain can come from 
archiving the hash values for each entity. 

On a technical note, we find a slight speed-up associated with implementing 
the field as a two-level dictionary. The upper level keys are tuples of Lorentz 
indices. These each index a lower level dictionary of all entities with the ap- 
propriate Lorentz structure. 

If, on inserting a new entity into a field structure, an equivalent entity exists, 
rather than adding the new item to the field its amplitude is merely combined 
with that of the existing entity. If the new amplitude + = 0, the entry 
is removed. This test is robust for integer arithmetic. Otherwise, the absolute 
value is compared to some tolerance, e.g. 10~^. We may worry about rounding 
errors: near fcj = the reduced vertex monomials are all adding in phase, and 
the deleted amplitudes may add to give a significant contribution. We therefore 
use a tolerance smaller than is finally required for the path expansion. As a 
last step only, we apply the final tolerance cut. The numerical values of the 
intermediate and final tolerances are found by trial and error, looking for 
robustness in the number of terms in the expansion. It is worth pointing out 
that Python carries out all floating point arithmetic in double precision. 

In the gluonic case of closed, traced contours the endpoint y is physically 
irrelevant. By ignoring it, we can identify more entities that are equivalent and 
further reduce the dictionary size. To find these, for each entity we impose x — 
y = whilst translating all the v's such that Vi = e^^. The field dictionary 
is then rehashed, looking for newly equivalent terms. Of course, care must be 
taken only to do this as a final step, and not to then multiply such contours 
together. 

As discussed in Section 2.5, a recursive action definition is more compact 
and less error-prone. Each smearing definition (such as "APE" or "ASQ") is 
predefined, and labelled by a unique string. Each path in the action is specified 
by three items: the coupling for this path, a list of signed directions and a list 
of link improvement method names. A full action (gluonic or fermionic) is 
specified by a list of such path specifications. For each path, the expansion 
routine is executed recursively to implement Eq. (27), converting it into a field 
that is fed to the next level of the nested link improvement. Finally, these fields 
can be manipulated or combined, before being output. 
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In some cases we need to combine complex conjugate monomials as per Eq. (21). 
This is most easily done by noting that E = 2ii^i.eai ~ E*. We loop over all 
unprocessed entities E E F, inserting —E* into F and marking E as now 
processed (and doubling its amplitude). This marking can be done by adding 
a Boolean flag to the entity data structure in Appendix C. 

Once we have the monomials, we have a choice of output format dictated by 
whether the output is to be read into the Feynman diagram evaluation code at 
compile time or run time. The former has the advantage that the compiler may 
be able to optimise the construction of the vertex functions; the disadvantage 
is that the size of the vertex functions may lead to code that is too long for 
the compiler to handle. Reading in the data at run time avoids this, but may 
in principle lead to slower code. In either case the Python code can be easily 
adapted to produce the correct output format. 

As an example, we describe the run time case. The output consists of a single 
ASCII file for each order of the perturbative expansion. Each file contains 
multiple entries, which could be a single line, or multiple lines for clarity. 
The entry contains the information in a single entity as whitespace separated 
values. For later storage in the Feynman diagram code it is convenient if 
the Lorentz indices are represented as a single integer in base D: n{fi) = 
J2l=i{lJ-i — 1)-D*~^. It is also useful if the entries for given /i arc consecutively 
numbered (although the order does not matter). The file is terminated by a 
blank entry with n(/i) = —1. 

In the Feynman diagram code a set of arrays should be defined to hold the 
vertex function data. For languages without allocatable arrays, we can arrange 
for the Python to write a set of compile time header flies that create arrays 
of the correct dimensions for a given set of vertex functions. 

We use a set of Fortran90 modules to read in the data flies, which can serve 
as a template for other languages. 



5 Conclusions 

Simulation of Symanzik and radiatively improved lattice field theory actions 
has become very popular in recent years. Associated renormalisation factors 
(and, indeed, the radiative improvement itself) can be systematically calcu- 
lated using lattice perturbation theory. The complicated nature of the im- 
proved actions and operators has, however, contributed to a backlog in this 
perturbative renormalisation programme. 

Having a flexible method for generating the Feynman rules automatically is 
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crucial to overcoming this backlog, and permitting a greater range of renor- 
malisation factors to be calculated. This paper provides just such a method, 
that is well suited to expanding all sectors of lattice QCD: gluons, relativistic 
fermions and heavy quarks. In addition to the Taylor expansion algorithm, an 
efficient implementation in the Python programming language is described, 
exploiting useful features of this language. 

Particular strengths of this algorithm include coping with arbitrary spin and 
colour trace structures in the action, allowing a nested definition of link im- 
provement and an intuitive way of defining the action to be expanded. 

The code is also very fiexible, and can be adapted to deal with most wrin- 
kles met in perturbative expansion. The first of these is tadpole improvement. 
Tadpole (or mean field) improvement aims to speed the convergence of per- 
turbation theory through dividing each link in the action by a factor u [27]. 
We can use perturbation theory to calculate u — 1 + diag + ■ ■ ■ as an expan- 
sion in the coupling, and treat the quantum effects as radiative counterterms 
in the action with couplings ritdi^asY (plus combinatoric factors), where rit 
was the number of factors of links in the path. We can expand the action as 
before, but must now do a separate expansion of the radiative counterterms. 
Consider building an action from links smeared as in Eq. (26). When we ex- 
pand W^{x)Wfj,{x + e^) we will get a link combination [z/, /x, — i/, i/, i/], and 
the question is whether such a term should be given a tadpole improvement 
factor of wr^ or -u^®. The former is more in keeping with the philosophy of 
mean field improvement, but in a simulation the latter is more convenient. 
The procedure, of course, is that the perturbative action should follow what 
is used in simulation. Either convention can be followed by assigning to each 
entity knowledge of the full length of the path from which it is derived. In the 
latter case of no cancellation, these lengths simply add on entity multiplica- 
tion. In the former case, some directional knowledge must be maintained to 
allow factors of u to cancel. Terms with different numbers of tadpole factors 
should be grouped separately; for this reason rit should be included in the key 
in the Python implementation of the entities and fields. 

This adaptability also makes the expansion algorithm described here useful in, 
for instance, chiral perturbation theory [28] or the double expansion needed 
for stochastic perturbation theory [6,29]. 

By describing and making available these algorithms and tools, we hope that 
lattice field theory calculations can reach a point where the choice of lattice 
action is not constrained by the availability (or not) of renormalisation con- 
stants. 
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A Spin algebra 

Each Wilson hne in an action has an associated spin matrix. Where this is 
not uniformly unity, wc must keep track of which spin matrix applies to each 
term in the reduced vertex. We do this by adding a single integer label to the 
entity list, r. For Dirac gamma matrices < r < 15, whilst for Pauli sigma 
matrices < r < 3. By convention, r = is the identity. 

When two spin factors are multiplied, the product is proportional to a single 
element of the spin algebra: 



(with no sum over k). This reduction can be easily encoded in the PYTHON 
vertex generation code through a small dictionary where each key (r^r-') 
indexes a list: [r'^,£fc]. 



B Pattern lists 

A Wilson line may be composed of a number of parts to which separate 
colour traces may have been applied. This will affect the value of the as- 
sociated Clebsch-Gordan coefficient. Physical examples include the traceless 
field strength operator, U^i, — Trf/^j,, and the adjoint Yang-Mills action, 
Tt Ufj^iyTi U^j^. In the former case, for instance, second order monomials either 
have colour structure of the form (TaT(,)cd or Tr(TaT;,) 5cd = — | Sab ^cd- 

As we do not calculate the Clebsch-Gordan coefficients in the vertex generation 
program, we need a method for distinguishing whether given gauge potentials 
within an entity are untraced, all traced together or in separate colour traces. 
This distinction is also important in ensuring we do not compress together 
entities with different colour structures. We distinguish these cases by adding 
an extra entry to the entity called a pattern list. 

A pattern list of order r is a; = (a;i, . . . ,0;^). Each positive integer element 
uji is associated with the corresponding factor of the gauge potential Aj. If Ai 




(A.l) 
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has not been traced, Ui — 0. All gauge potentials with the same value of the 
pattern component are understood to be contained in a single colour trace. 
For instance, AiA^ Tr(y42v44) would have a pattern list (0, 1, 0, 1). 

Applying a colour trace to an entity modifies only the pattern list u? — > w', 
with 

I l + max(a;i,...,a;^) , = , 
= \ (B.l) 

We stress that the actual value of cUj has no meaning. It is therefore convenient 
to arrange at all stages that the first non-zero element in the list is 1, the second 
2 etc. Taking the example above, TrfAiAg Tr(A2A4)] = TT{AiAs)Tr{A2A4) 
has a relabelled pattern list (1, 2, 1, 2). Gauge invariance precludes application 
of a trace to entities for which x ^ y (i.e. to contours which are not closed). 
The SU{N) generators are traceless, so when a colour trace applies to only 
one factor of the gauge potential, we may delete the entity (such as when 
tracing an entity of order r = 1). 

We define multiplication u?^ * u?^ by: 

{ul . . . ,<) * {ul . . . = (^, . . . . . . ,<) (B.2) 

where 

[ ui + max(a;i, . . . , ul.) , ui . 

When all elements of w are not 1, the symmetries of Eq. (23) are not present. 
The group of symmetrisation operations carried out in the vertex generation 

code must therefore be reduced from Zr, or it may be simpler to postpone 
all symmetrisation until the Feynman diagram are evaluated for specific mo- 
menta. 

We note that pattern lists can also be used to label the taking of real or 
imaginary parts in an action in a similar way, using positive and negative 
entries to distinguish the two. 



C Entity algebra 



Taking into account Appendices A and B, an entity E consists of 

E^{fx;x,y;vi,...,Vr;oj;T;f). (C.l) 

The colour trace pattern list uj and spin index r are optional and need not be 
included in all situations. The start site x may also be omitted (see Section 4). 
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The complex conjugate entity E* has amphtude (—1)''/* and the sign of all 
D-vectors reversed. 

Multiplication by a scalar p G C acts only on the amplitude: 

pE = {fi; X, y;vi,..., Vr] uj; r; pf) . (C.2) 

We translate an entity by D-vector c G A using: 

T^E = (ni, ...,Hr;x + c,y + c;vi + c,...,Vr + c;u;; r; /) (C.3) 

Two entities are said to be equivalent if the lists can be rendered identical 
under a translation and rescaling: 

E' = E^ iff 3 c G A, p G C s.t. T^E' = pE^ . (C.4) 

Non-commutative multiplication of two entities is defined by: 

E' = EU E^ = {fi'; x\ + C; 

...,<, + C, ... , vi^ +C;cvU lj'; t'; f) , (C.5) 

i. e. the path from which the second entity was derived is first translated by a 
vector C = y^ — x\ to start where the first finished. The resulting entity is 
of order r = Vi + rj. The Lorentz list /-i' is the concatenation of lists /j,^ + /j,^ . 
The spin indices yield rV-' = efcx'^ as per Eq. (A.l). Note the amplitude 
/' = Ek ^Cr^ /* contains a combinatoric factor arising from having separated 
out the (r!) Taylor expansion factors from the amplitude in Eq. (8). 

The action of the permutation operator a = {ai, a2, ■ ■ ■ , Cr) on a list I yields 
(/o-j, . . . , Ifj^). We can apply it to entities which are closed x = y, simply traced 
{uji = 1 V i) and where the real part has been taken (to ensure Eq. (22) holds) 

a ■ E = {a ■ fi;x,x;v^^,. . .,v^/,a ■ u;;t; xM)!) , (C.6) 
noting in this case that a ■ uj = uj. Eq. (23) defines Xr{^)- 

By extension Eqs. (C.2, C.3, C.6) and the colour trace are applied to a field 
F by operating on each of its constituent entities E & F. 
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